
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the field caused by point charges.
%   xx_trg   = locations of target points.
%   xx_src   = locations of sources points.
%   qq       = the strength of the point charges.
%   kh       = helmholtz parameter.
%
% The output u is the vector
%
%   u(i) = sum_{j} qq(j) * H_0(kh*|xx_trg(:,i) - xx_src(:,j)|)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u = pot_from_charges(xx_trg,xx_src,qq,kh)

% number of points on the target contour
ntrg  = size(xx_trg,2);
% number of source points
nsrc  = size(xx_src,2);
% dist1 = x_trg(i) - x_src(j)
dist1 = (xx_trg(1,:)' * ones(1,nsrc)) - ones(ntrg,1)*xx_src(1,:);
% dist2 = y_trg(i) - y_src(j)
dist2 = (xx_trg(2,:)' * ones(1,nsrc)) - ones(ntrg,1)*xx_src(2,:);
% dist =  | x_trg(i) - x_src(j) |
dist  = sqrt(dist1.*dist1 + dist2.*dist2);

% u(i) = sum_j H1(h | x_trg(i) - x_src(j) |)* qq(j)
u = besselh(0,kh*dist)*reshape(qq,numel(qq),1);

return
